In [1]:
from __future__ import print_function

In [2]:
from functools import partial

%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import astropy.units as u
import numpy as np

from astropy.modeling import models, fitting

Now we try and build the mfe single tube model.


In [3]:
import pysac.mhs_atmosphere as atm
from pysac.mhs_atmosphere.parameters.model_pars import mfe_setup as model_pars


WARNING:traits.has_traits:DEPRECATED: traits.has_traits.wrapped_class, 'the 'implements' class advisor has been deprecated. Use the 'provides' class decorator.

In [4]:
# Cheeky Reset to Photosphere
#model_pars['xyz'][4] = 36641.22*u.m
#model_pars['xyz'][5] = 1587786.3*u.m
#model_pars['B_corona'] = 0.00055*u.T
#model_pars['radial_scale'] = 0.03938*u.Mm
##model_pars['radial_scale'] = 0.05*u.Mm
model_pars['pBplus'] = 1.2e-3*u.T
#model_pars['chrom_scale'] = 0.4*u.Mm

In [5]:
keys = ['pBplus', 'B_corona', 'coratio', 'chratio', 'radial_scale', 'chrom_scale', 'corona_scale', 'xyz']
for key in keys:
    print(key, model_pars[key])


pBplus 0.0012 T
B_corona 0.00055 T
coratio 0.0
chratio 1.0
radial_scale 0.03938 Mm
chrom_scale 0.4 Mm
corona_scale 0.25 Mm
xyz [<Quantity -1.0 Mm>, <Quantity 1.0 Mm>, <Quantity -1.0 Mm>, <Quantity 1.0 Mm>, <Quantity 0.036641221 Mm>, <Quantity 1.5877863 Mm>]

In [6]:
# model setup
scales, physical_constants = atm.units_const.get_parameters()
option_pars = atm.options.set_options(model_pars, False, l_gdf=True)
coords = atm.model_pars.get_coords(model_pars['Nxyz'], u.Quantity(model_pars['xyz']))

#interpolate the hs 1D profiles from empirical data source[s]
empirical_data = atm.hs_atmosphere.read_VAL3c_MTW(mu=physical_constants['mu'])
table = atm.hs_atmosphere.interpolate_atmosphere(empirical_data, coords['Zext'])


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-e26610594ef7> in <module>()
      2 scales, physical_constants = atm.units_const.get_parameters()
      3 option_pars = atm.options.set_options(model_pars, False, l_gdf=True)
----> 4 coords = atm.model_pars.get_coords(model_pars['Nxyz'], u.Quantity(model_pars['xyz']))
      5 
      6 #interpolate the hs 1D profiles from empirical data source[s]

/opt/miniconda/envs/thesis/lib/python2.7/site-packages/pysac/mhs_atmosphere/parameters/model_pars.pyc in get_coords(Nxyz, xyz)
    144     dz=(xyz[5]-xyz[4])/(Nxyz[2]-1)
    145     Z    = np.linspace(xyz[4],
--> 146                         xyz[5],Nxyz[2])
    147     Zext = np.linspace(Z.min()-4.*dz, Z.max()+4.*dz, Nxyz[2]+8)
    148     coords = {

/opt/miniconda/envs/thesis/lib/python2.7/site-packages/numpy/core/function_base.pyc in linspace(start, stop, num, endpoint, retstep, dtype)
    113 
    114     if endpoint and num > 1:
--> 115         y[-1] = stop
    116 
    117     if retstep:

ValueError: setting an array element with a sequence.

In [ ]:
table['rho'] += table['rho'].min()*3.6

In [ ]:
#==============================================================================
#calculate 1d hydrostatic balance from empirical density profile
#==============================================================================
# the hs pressure balance is enhanced by pressure equivalent to the
# residual mean coronal magnetic pressure contribution once the magnetic
# field has been applied
magp_meanz = np.ones(len(coords['Z'])) * u.one
magp_meanz *= model_pars['pBplus']**2/(2*physical_constants['mu0'])

pressure_z, rho_z, Rgas_z = atm.hs_atmosphere.vertical_profile(coords['Z'], table, magp_meanz,
                                                               physical_constants, coords['dz'])

In [ ]:
x, y, z = u.Quantity(np.mgrid[coords['xmin']:coords['xmax']:1j*model_pars['Nxyz'][0],
                              coords['ymin']:coords['ymax']:1j*model_pars['Nxyz'][1],
                              coords['zmin']:coords['zmax']:1j*model_pars['Nxyz'][2]], unit=coords['xmin'].unit)

In [ ]:
xi, yi, Si = atm.flux_tubes.get_flux_tubes(model_pars, coords, option_pars)

In [ ]:
Si[0][0] = 0.1436*u.T

In [ ]:
allmag = atm.flux_tubes.construct_magnetic_field(x, y, z,
                                                 xi[0], yi[0], Si[0],
                                                 model_pars, option_pars,
                                                 physical_constants,
                                                 scales)
pressure_m, rho_m, Bx, By ,Bz, Btensx, Btensy = allmag

In [ ]:
# local proc 3D mhs arrays
pressure, rho = atm.mhs_3D.mhs_3D_profile(z,
                                          pressure_z,
                                          rho_z,
                                          pressure_m,
                                          rho_m)
magp = (Bx**2 + By**2 + Bz**2) / (2.*physical_constants['mu0'])
energy = atm.mhs_3D.get_internal_energy(pressure, magp, physical_constants)

In [ ]:
x_lin = np.linspace(coords['xmin'], coords['xmax'], model_pars['Nxyz'][0])
bmag = np.sqrt((Bx**2 + By**2 + Bz**2))[:,64,0].to(u.mT)

gaussian = models.Gaussian1D()
bmag_fit = fitting.LevMarLSQFitter()(gaussian, x_lin, bmag)

fwhm = 2.*np.sqrt(2*np.log(2))*bmag_fit.stddev.value
fwhm = np.abs(fwhm) * u.Mm

fwtm = 2.*np.sqrt(2*np.log(10))*bmag_fit.stddev.value
fwtm = np.abs(fwtm) * u.Mm

u.Quantity([fwhm, fwtm]).to(u.km)

In [ ]:
%matplotlib inline
fig = plt.figure(figsize=(15,8))
plt.plot(x_lin, bmag.to(u.mT), 'x')
X = np.linspace(coords['xmin'], coords['xmax'], 1000)
plt.plot(X, bmag_fit(X.value))

plt.xlabel("X [{}]".format(x.unit))
plt.ylabel("Magnetic Field Strength [{}]".format(bmag.unit))
plt.axhline(y=bmag_fit.amplitude.value/10)

Now for some plotting


In [ ]:
import yt

In [ ]:
def magnetic_field_strength(field, data):
        return np.sqrt(data["magnetic_field_x"]**2 + data["magnetic_field_y"]**2 + data["magnetic_field_z"]**2)
yt.add_field(("gas","magnetic_field_strength"),
             function=magnetic_field_strength,
             units = "T")

In [ ]:
def plasma_beta(field, data):
        return data['pressure'] / data['magnetic_pressure']
yt.add_field(("gas","plasma_beta"),
             function=plasma_beta,
             units = "")

In [ ]:
def alfven_speed(field, data):
    return np.sqrt(2.*data['magnetic_pressure']/data['density'])
yt.add_field(("gas","alfven_speed"),
             function=alfven_speed,
             units = "m/s")

In [ ]:
def sound_speed(field, data):
    cspeed = np.sqrt(physical_constants['gamma']*pressure/rho)
    return cspeed
yt.add_field(("gas","sound_speed"),
             function=sound_speed,
             units = "m/s", force_override=True)

In [ ]:
bbox = u.Quantity([u.Quantity([coords['xmin'], coords['xmax']]),
                   u.Quantity([coords['ymin'], coords['ymax']]),
                   u.Quantity([coords['zmin'], coords['zmax']])]).to(u.m).value

In [ ]:
data = {'magnetic_field_x':yt.YTQuantity.from_astropy(Bx.decompose()),
        'magnetic_field_y':yt.YTQuantity.from_astropy(By.decompose()),
        'magnetic_field_z':yt.YTQuantity.from_astropy(Bz.decompose()),
        'pressure': yt.YTQuantity.from_astropy(pressure.decompose()),
        'magnetic_pressure': yt.YTQuantity.from_astropy(magp.decompose()),
        'density': yt.YTQuantity.from_astropy(rho.decompose())}

ds = yt.load_uniform_grid(data, x.shape,
                          length_unit='m', magnetic_unit='T', mass_unit='kg',
                          periodicity=[False]*3, bbox=bbox)

In [ ]:
ds.index.grids[0]['density'].min(), ds.index.grids[0]['magnetic_field_z'].max(),

In [ ]:
ds.index.grids[0]['pressure'].min()

In [ ]:
ds.index.grids[0]['magnetic_field_z'][64,64,:]

In [ ]:
slc = yt.SlicePlot(ds, 'x', 'alfven_speed', origin='lower-center-domain', axes_unit='Mm')
slc.set_cmap('all', 'OrRd')

seed_points = np.zeros([11,2]) + 1.52
seed_points[:,0] = np.linspace(-0.99, 0.95, seed_points.shape[0], endpoint=True)
slc.annotate_streamlines('magnetic_field_y', 'magnetic_field_z', field_color='magnetic_field_strength',
                         plot_args={'start_points':seed_points, 'density':15, 'cmap':'winter',
                                    'norm': mpl.colors.LogNorm(*ds.all_data().quantities.extrema("magnetic_field_strength"))})
slc.annotate_contour('plasma_beta', clim=(1,50), label=True, take_log=False)

In [ ]:
slc = yt.SlicePlot(ds, 'x', ['pressure','magnetic_field_strength'], origin='lower-center-domain',
                   axes_unit='Mm')
slc.set_cmap('all', 'BuPu')
#slc.annotate_contour('plasma_beta')

seed_points = np.zeros([11,2]) + 1.52
seed_points[:,0] = np.linspace(-0.99, 0.95, seed_points.shape[0], endpoint=True)
slc.annotate_streamlines('magnetic_field_y', 'magnetic_field_z', field_color='magnetic_field_strength',
                         plot_args={'start_points':seed_points, 'density':15, 'cmap':'winter', 'linewidth': 2,
                                    'norm': slc.plots['magnetic_field_strength'].image.norm})
slc.annotate_contour('plasma_beta', label=True, take_log=False, clim=(1,50))

In [ ]:
slc = yt.SlicePlot(ds, 'z', 'magnetic_field_strength', origin='center-domain',
                   axes_unit='km', center=("max", "magnetic_field_strength"), width=(0.5, "Mm"))
slc.set_cmap('all', 'Oranges')
slc.annotate_contour('magnetic_field_strength', take_log=False,
                     plot_args={'levels':[ds.all_data().quantities.extrema("magnetic_field_strength")[1]/2.],
                                'linewidths':3,
                                'colors':'black'})

In [ ]:


In [ ]: